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ABSTRACT 

The application of Boundary integral Equations to elastoplastlc problems is re- 
viewed. Details o r the analysis us applied to torsion problems and to plane problems 
is discussed. Results are presented for the elastoplastlc torsion of a square cross 
section bar and for the plane problem of notched beams. A comparison of different 
formulations as well as comparisons with experimental results are presented. 

INTRODUCTION 

Methods of analysis in elasticity and plasticity, ns in most other scientific and 
engineering fields, have been revolutionized by the advent of the modern digital com- 
puter. Thus tiie nvailability of the computer made it possible to implement practically 
purely numerical methods such as finite difference and finite elements, as well as 
analytical methods such as complex variable methods. 

Similarly the boundary integral equation (BIE) methods, while having their origin 
in classical elasticity, have only in recent years begun to play a significant role in 
solid mechanics. Solutions to problems in elasticity by these BIE methods have been 
obtained by various investigators using different formulations, ns for example in Refs. 
(J0). A review o' much of the literature is given in Ref. 0. 

The extension of the BIE method to elastoplastic problems has received much less 
attention, The basic theories and equations have been ‘cumulated in Refs. 0.4) but 
few applications have been reported. The present paper reviews some of these appli- 
cations and presents details of the analyses as applied to elastoplastic torsion prot>* 
lems, and to the plane elastoplnstic problem, with particular reference to edge-notched 
beams in bending. Comparisons of different formulations as well as comparisons with 
experimental results are presented. 


ELA8TOPLA8TIC TORSION 


The clastcplaetic tom Ion problem can be formulated in several ways (Ref. 0). For 
example, using the Prandtl stress function, F, the basic differential equation becomes 


V 2 F = -2Go - 2G 



s «*,y) 


(1) 


where G is the shear modulus, o^ the linear coefficient of thermal expansion and 
4 ax’* sy Bre pl®*Uc ®hear strains. The corresponding B1E is 


irF(p) = ff f(<J)ln 
R 



F’ In r. 




(2) 


where primes denote derivatives with respect to the outward normal and where the 
coordinate system and the associated notation are shown in Figs. 1 and 2. For an in- 
terior point P, the multiplier of F(p) in Eq, (2) becomes 2ir instead of ir. 

As an illustration of the use of Eq. (2) and its ability to solve the elastoplastic 
torsion problem, consider the case of a circular shaft of radius a. The radial coor- 
dinate will be designated by p, to distinguish It from r, the distance between the 
fixed point and the variable point appearing in Eq. (2). In polar coordinates, because 
of symmetry, the function f appearing in Eq. (2) becomes 


f = -2Go 


+ 2G JJ. / p \ 


For linear strain hardening 


c z o = Ap ^ B 


( 3 ) 
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where 


A = 


2 


3o 


m 

3 + 2(1 + v) 

t - m 
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3 + 2(1 + i>) — 22— 
1 - m 


v is Poisson's ratio, c Q Is the yield strain, and m is the strain hardening parameter 
(slope of strain hardening curve divided by the modulus). 

On the boundary F(a) = 0 and because of axial symmetry F'(a) = constant. Eq. (2) 
then becomes 


0 “ 2G 


f/t - a + — ^ln dx dy - F'(a) j In r N dq 


which ifion solving for F'(a) gives 

F'(a) - G|a<2A - o) +2Bj 


(5) 


«J> 


.Hence for any interior point 


2irF(p) - 2G /ft-'-) mr pQ dxdy. V'wf c 

,n r (jq tIq 


or 


K(p) ~ | |(2A - o)(p 2 - a 2 ) + 4D(p - a)J 
and the shear stress t is given by 

r -?- 2c [(f • a ) ,p - b; 1 


(7) 
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which agrees with the solution obtained in a entirely different fashion in Ref. (h). Note 
that this solution is valid only in the plastic region, that is, for p <p c , where p^, 
the elastic plastic boundary, is given by 


P 
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2(1 + r)e o 


(9) 


The formulation given by Eq. (2) can of course be used to obtain the elastoplastic 
solution for almost any shape cross section and any type of strain hardening. In gen- 
eral, however, it would seem that a tormulation in terms of the warping function (axial 
displacement) should be preferable, since the warping function is physically more 
meaningful than the stress function and more importantly, the distinction between sim- 
ply connected and multiply connected regions disappears. We will therefore formulate 
the problem in terms of the warping function and show in some detail how the solution 
is obtained for a bar with a square cross section. 

The BIE in terms of the warping function w is (Ref. (0)) 


where 


irw(p) * j f f(Q)l» r pQ dA t ^w<q)(ln r^y dq - w'(q)ln r^ dq 
R q 
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(11) 


Note that Eqs. (U and (10) imply the use of the deformation theory of plasticity. How- 
ever, as .diawn by Prager (Ref. (7)(, both the total and incremental theories of plas- 
ticity furnish the same solution to the torsion problem provided either the cross section 
Is circular or the material is perfectly plastic. It is reasonable to assume, therefore, 


(hat this will be approximately true (or most practical problems. Indeed, it has been 
shown in Ref. 0 that for the case of a square cross section with strain hardening 
there is little difference between incremental and deformation theories, the use of 
inc remental theory does not appreciably complicate the problem and can be used in a 
stepwise manner as for the plane problems to be discussed subsequently. 

The boundary condition lo be satisfied ly the warping function is (Ref. 0) 

— ® w* * o(fy - mx| + 2 (* l xz + mt vr.) (* 2 ) 

where £,m are die direction cosines of ihe outward normal, and where for a rectangular 
boundary segment parallel to one of the coordinate axes, the second term on the right 
side of Eq (12) always vanishes. The normal derivative of the warping function appear- 
ing in Eq. (JO) is ihus known from Eq- (12) and the only unknown In Eq. (10) is w(p). 

Numerical Procedure 


To solve Eq. (10) for the unknown function w(p), the straightforward procedure of 
replacing the integrals by summations cun be used. The boundary is divided into n 
intervals with a nodal point taken at the center of each interval. The unknown function 
is ussumed constant over cuch Interval, Similarly, the region R is divided Into a 
number of colls and the function f assumed constant over each cell. Eq. (10) is then 
written for each nodel point as follows: 


n n 

EvVvEV! 

FI FI 


R. 


t ^ 1,2 n 


(13) 


The coefficients tty, by, and Hj ttre given In Appendix A. We thus have n equations 
for the n unknowns, w^. This set of equations can readily be solved by any standard 
procedure. 

Once the Wj are known on the boundary, Eq. (10) ean be used to calculate w at 
any interior point, with it replaced by 2 it. However, in order to calculate the strains 
the derivatives of w are needed. These can be obtained by differentiating Eq. (10) 
under the integral sign lo give 
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For 3w<P>/0y wo Interchange x and y. 

Again the integrals are replaced by sums resulting in 


x 2x lfe 


<Vk A uk + w k m k B ijk“ w fc c iik> + 


^ f W D iJklj 


where the coefficients Ay^, By^, Cy k , and Dy^j are listed in Appendix A. £ ie 

k| i 

the sum for all the plastic cells in the region. 

From the derivatives of w the total strains are computed as 


= 1 Ly *SSL\ 
2 \ ax) 

,.iU+sx) 

' 3 \ By) 


The plastic strains appearing in the definition of the function f(x,y) are of course 
in turn nonlinear functions of the warping function w. They can bo determined from 


e P =_H f 
xz 4 xz 
e et 


f P =— E. c 


where 


eet ^V (cx/+<£ - )2 


( p- f(t ct> 

Eq. (19) represents the uniaxial stress-strain curve in terms of equivalent plastic 
strain against equivalent total strain; that is, 

1 a e 

where <r Q , the equivalent stress, represents the stress on the uniaxial stress-strain 
curve and e the plastic strain on tnat curve. Thus, for a given stress-strain curve, 
the relation between Cp and represented by Eq. (19) can be determined using 
Eq. (20). For the case ol linear strain hardening, the relation (19) can be written as 


‘cf^ 1 * 

P 1 + - (1 + i') — — 
3 1 - m 


( 21 ) 


Since as noted the (unction f is implicitly a nonlinear function of the warping func- 
tion w, the solution is obtained iteratively by starting with f * 0 and calculating im- 
proved values via Eqs. Q3), (is), (16), (17), (il), and back to (13), etc. This Is the 
method of successive elastic solutions or method of initial strains. The tangent modu- 
lus method could be used equally well and may save on computer time. Complete do- 
tails of the calculations performed herein are given in Ref, (6). 

Results 

Calculations were performed by this technique for a bar of square cross section as 
shown in Fig. 3. The dimensionless angle of twist per unit length i), defined as oa/< Q ; 
where a Is l/2 the side of the square, was increased in steps of one from ) * 1 to 
1 = 6. Linear strain hardening was assumed with values of the strain hardening param- 
eter taken as 0 (perfect plasticity), 0.05, 0.1, and 0.2. Poisson's ratio waB assumed 
as 0. 3 in all calculations. 

For 3 * 1, the bar is clastic and a comparison was made between the analytical 
solution as given, for example, in Ref. (g) as well as with the finite difference solution 
of Ref. (10). The results are shown in Tables 1 to III. Tabic 1 shows the warping func- 
tion as computed on the boundary of the bar cross section. The comparison with the 
analytical solution of Ref. 0 shows very good agreement with Just four unknowns to 
solve for in the boundary integral method. 

Table II show t the comparison for the maximum shear stress (at the center of the 
edge of the square) and the moment with the analytical solution of Ref. (0) n^d the finite 
difference solution of Ref. (JW). Again it is seen that with just four unknowns in the 
boundary integral method very good results arc obtained, as good as the results ob- 
tained for the finite difference method using 55 unknowns. 

Table III presents the dimensionless shear stress distribution in the x-direction 
(T J<s ,/2Ge 0 ) throughout the cross section using 10 unknowns for the boundary integral 
method and 55 unknowns for the finite difference method. Again excellent agreement 
was obtained. Actually, the results with four unknowns using the boundary integral 
method are almost as good, but the results with 10 unknowns are presented to match the 
actual (x,y) values of the finite difference results without having to cross plot. 

The dimensionless angle of twist per unit length il was then increased in unit stopB 
to a maximum value of d - 0 for each value of tire strain hardening parameter m. The 
total boundary was divided into 80 intervals resulting in 10 equations for 10 unknowns. 
Several test calculations were made with fewer intervals, and the results indicated that 
using 48 intervals (six unknowns) changed tire moment and maximum stress by at most 
one in the third significant figure and changed the maximum plastic strain by nbout 3 per- 
cent. All the subsequent results are therefore shown for 80 intervals (10 unknowns), 
although from an engineering viewpoint 48 or even 32 intervals would be sufficient. 

The results of the calculations are summarized in Tabic IV and Figs. 1 to 6. 

Fig. 4 shows the dimensionless moment defined as M* - M/2Gi Q a' for various values 
of 3 and m. Fig. 5 shows the corresponding dimensionless maximum shear stresses 
defined as r max = r/2Ck 0 and Fig. 0 show's the spread of the plastic zones with an in- 
crease of tile angle of twist J. 

The degree of convergence of die iterative process was determined from a relation 
of the form 


( 32 ) 


i=i 

where M le the total number of point* P ( (lowing plastically and k-1 and k arc two 
successive iterations. The convergence criterion K can be made as small as desired. 
In all the calculations (he convergence number K was taken cither as 0. 0001 or 
0.00001. In many of tho calculations both numbers wero used in turn. The differences 
in the results were found to be insignificant. For example, the number of iterations 
for convergence for the case of maximum plastic flow, which occurred for j = fi and 
m = 0, was 39 for K = 0. 0001 and S3 for K = 0.00001, and the results were all the 
same to at least throe significant figures. For the case 3-5 and M = 0, the number 
of iterations for K = 0 0001 was 33. For the same case using finite differences, 

203 iterations were required. 

The boundary integral method is thus seen to be very suitable for the elastoplastic 
analysis of the torsion of prismatic bars. Very good accuracy can be obtained by using 
relatively small seta of linear algebraic equations. 

A comparison with the finite difference method indicates a groat savings in the 
number of unknowns that have to be determined and also a much faster convergence rate 
using the method of successive elastic solutions for both formulations. This should be 
reflected in appreciable savings in computer time, although computer time is not a 
limiting factor in any case for the torsion problem. 

The boundary Integra- ncthod can readily be programmed in a straightforward 
manner for a digital computer. The use of the warping function to formulate the prob- 
lem permits applying the method with equal ease to both simply connected and multiply 
connected bodies. 

THE PLANE PROBLEM 

As for the torsion problem, the plane problem can be formulated In several ways, 
as a nonhomogeneous biharmonic problem for the stress function (Ref. (3)), or in 
terms of the Navier equations of equilibrium for the displacements (Refs. 5, 4), Both 
methods will be applied herein to the problem of an edge-notched beam in pure bending. 


The Biharmonic Formulation 


The problem of determining the state of stress and strain Ir. a plane elastoplastic 
problem cun be reduced to solving the following inhomogcncou' biharmonic equation for 
the Airy stress function, v>, as shown in Ref. (5) 

V*c> = g(x, y) (23) 


g(x, y> - - 


K 

" 2 

, 2 

2 

1 ~jj 

L>y 




( 24 ) 




tor the piano strain case, and 


tf<x,y) 



( 25 ) 


for the plane tress ease, where t^, and tjjy rqjrosont the accumulation of plas- 
tic strain increments from the buj'i.tning of die loading history up to, but not Including 
the current increment of the load, and Ac^, AcJ, and Ac^ are the increments of 
plastic strain due to the current increment of load. 

The stress function >/> must satisfy appropriate boundary conditions. Tor the 
problem under consideration (Fig. 7), i?(x,y) and its outward norma? derivative dip/Qa 
must satisfy the following boundary conditions (Ref. (11) : 


ip(x, v) = 0, — = 0 along boundary OA and OA' 
un 

^(x,y) = 0j = 0 along boundary AB and A'B' 

an 


^<x,y) - 4 7) 4 ° m ax(j 4 “ x 4 7) » 


— - 0 along boundary BC and B'C 


yn 




2 


^(x,y) - mitx — ; ~ = 0 along boundary CD and CD' 
tt yn 


(26) 


To solve Eq. (23) by means of the boundary integral method, use is made of 
Green's second theorem to reduce this equation to coupled integral equations, as shown 
in Refs. (3) and (11). The result is 


8iW>(x,y) - ff pgU.n)d4 dij = / L>— ( v 2 (J) _ 

R Jc *- 8 " 


4ir</>(X, 


and 


H£v 2 ,j p 

dq 

for 

PC R 

yn an yn 
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il£ v 2 !) 1 + 

dq 

for 

PC C 

an an yn 



(28) 


2 ff*<x,y) - ff g( 4 , 7 j)ln r d£ eh) 

It 




(In r) - — in r 
yn 


dq for PC R (29) 


for PCC 


»*(x,y) 


f f gtf rOidn-l L — (In r) - ~ In r dq 

V Jr L a " im j 


< 30 ) 


where 


* * V 2 ^ 


p m r 2 In r 

and r(x,y; £ ,rj) is the distance between any two points P(x,y) and q(£,n) in (he re- 
gion R bourutod try the curve C, such that PC R + C and qc C (Fig. 8). 

Eq. (27) would, for a known function g(X,y), give us directly a solution to the 
bihnrmonlc Eq. (23) provided the functions v>(x,y), Ov><x,y)/0n, V 2 u5(x,y), and 
a[V 2 «5(x, y)]/i)n were known "a the boundary C. 

However, only the stress function v> and its outward normal derivative av>/on 
are specified (Eq. (26)). The values of v v>«* and a(V tP)/an*a*/an on the 
boundary must be compatible with the given values of and Ov>/gn. To assure this 
compatibility, we have to solve the system of coupled integral Eqs. (28) and (30), 
which contain the unknown functions ♦ and o4/un. 

Once the values of ♦ and a* /an on the boundary C of region R are known we 
can proceed with the calculation of the stress field in (he region R utilising Eq, (27) 
and tiic equations which define u>, namely. 


. uV r . _ 0 2 «g- 

ay 2 * ax 2 * 


(31) 


The calculation of the function g(x,y), which is obtained Iteratively, will be dis- 
cussed subsequently. 


Solution of Ihe integral Equations 

To solve the system of coupled integral equations analytically, a numerical method 
is utilized in which the integral Eqs. (28) and (30) arc replaced by a system of simul- 
taneous algebraic equations. 

For simplicity of notation the normal derivatives arc denoted by prime super- 
scripts. The boundary is divided Into n intervals, not necessarily equal, numbered 
consecutively in the direction ot increasing q. The center of each interval is desig- 
nated as a node. The values of * and are assumed constant on each Interval and 
equal to the values calculated at the node. 

In similar manner the interior of region R is covered by a grid, containing m 
cells. The cells do not have to have equal areas. Their nodel points are located at the 
centroids. The value of g(£,n) is assumed constant over cadi cell and equal to the 
value calculated at the centroid. The arrangement of boundary and interior subdivisions 
is shown in Figs. 9 and 10. 

Using these assumptions, Eqs. (28) and (30) can be replaced by a system of 2n 
simultaneous algebraic equations with 2n unknowns, that is, *j and +J. 
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(32) 


*v 
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m n 

In r |ktRA»K * E“u»i * vp 

k*l ]*1 


' J «“ (BA) k ;Z/Vj + Vi + V j + V1‘ 

k-l J*1 

where 1*1,2, 3,- . . , n, r [k in the distance from node to the eentrdd of the 
cell, is the a rat of the k dl cell, and 


<ln r ( j)* dq 


f, tJ ' jf 


b i j — lnr u dq 


C,j = J „{j dq 

d ij a - f »ll ikl 

e U - / < v2f> ij > * <kl 

f U ' * / * J 


(33) 


where integration is taken over the j dl interval, and r.^ is the distance from l 
node to any point in the J* interval. The normal derivatives in Eqs. (33) are taken 
on tile j dl interval. 

Tor curved boundaries Ihe coefficients given by Eqs. (33) can be evaluated, if nec- 
essary, by Simpson's rule for i 4 j. For i - j, because of the singular nature of the 
integrand, the integrals for the coefficients must be evaluated by it limiting process. 
For boundary intervals, such us lor the problem treated herein, which can bo repre- 
sented by straight lines a closed form solution cun be obtained for these coefficients. 

Boundary Eqs. (32) expressed in matrix form become 


i»ij- 6 ifi 

i-i.r 
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i ln r tkl 

loj |0J 

n x n 

n * n 


n > 1 


n • m 

n ■ n n * n 
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n » 1 


n ■ m 
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^ J 


- 
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m x 1 


V. 


U’jJ 

n x i 

l*jl 

n x 1 


l« 


(31) 


Thus, tho problem In reduced to the solution of tho following matrix system: 

[B]{Xj *{R} 


( 3 ® 


whore (Bj is 2n x an matrix and {X} and (R; are 2n x 1 column matrices. 

Matrix [Bj is dependent only on geometry, that is, number of nodes and their dis- 
tribution on the boundary. Since tho matrix {Rj contains the nonl near function g(£,rj), 
which depends on the stress field and therefore on matrix {X}, an iterative process will 
be used to obtain the solution. 

To calculate stresses, at any nodal point in the region R, from the stress function 
< 1 , wc need not perform any numerical differentiation. Eq. (37) can be differentiated 
under the integral sign and once * and ♦ ' arc known on the boundary the stresses 
can bo obtained hy the same typo of numerical Integration as In Eq*. (33). Applying 
Eqs. (31) to Eq. (27) yields for tiro case of a rectangular grid tho following stress 
equations: 



lL*zl£ + A 

(x-o 2 My ->)) 2 J, 


(ftA) 


k X <-Vj~ VJ f WVi' 


ik 

lA 




xy 


m 

(*,y), ^ [ ~ Al* {g A) k +■ £ <G ij v» j + IlyC) + I t j*j + Kjj*]) 

Z— fU^-4r + (y-n)J lk J=i 

iA 


where now 1 = 1, 2, 3,. . m refers to tire centroid of the i lil coll, 6 x and 6 
represent, respectively, x-direetional and y-direetional dimension of the coll. The 
coefficients Ay, By, Cy, By, Ey, Fy, Gy, 1 1 j j, ly, and Ky arc obtained by ap- 
propriate differentiation under the integral sign of the coefficients given by Eqs. (33) 
and are listed in Appendix B- 


The stress function <p Ik not constant on the loaded boundaries BC and B'C'. 

The assumption that U piece-wise constant may lead to appreciable errors In the 
numerical row u Its. To eliminate this source 1 of error, the summations given in 
Eqs. (.12) and (36) for intervaia tying on the loaded boundaries and involving the stress 
function are replaced by direct integration. 

Boundary Interval and Interior Grid Hlzs 

The number of nodal points prescribed for the boundary is theoretically unlimited 
However, computer storage capacity for the computer used and difficulties associated 
with Inversion of large matrices limited the order of die coefficient matrix [BJ of 
Eq. (35) used herein to HO, 

Because of geometric and loading symmetry about the x-uxis, it is possible to re- 
duce the number of unknowns. For 2n total number of nodul points the number of 
equations and unknowns <kj and $|, la reduced from 4n to 2n. Additional reduction 
in the number of unknowns is uecompllshcd by taking into consideration 81, Venant's 
effect at the loaded boundaries (Ref. (12)). 

Since the vicinity ol the crack tip Is of greatest interest, a fine nodal spacing along 
the notch was chosen. To reduce the error introduced by the change in the interval 
size (Ref. (13)) around boundary points A and A' and at die same time to obtain fine 
resolution at the tip of the notch, the boundary along the notch was divided into n number 
o' ,«!• ■ ils progressively decreasing in length. The rate of change in the interval 
f ngtl, and the resulting length of the smallest interval was feund to have a groat influ- 
ence on thv stress field in the vicinity of Uie tip of the notch. The rate of change in die 
interval's length along theso boundaries wus optimized by the method presented in 
Ref. (14) . For the cases considered optimum ratios of the lengdis of two consecutive 
boundary intervals were found to be in the range of 1.08 to 1.10. The resulting small- 
est dimensionless boundary interval length varied from 0.0001 to 0.0002. A set of 
140 equations containing 140 unknowns was used. Note 1 that the corner points are al- 
ways designated as interval points, never as nodal points, thus eliminating discon- 
tinuous functions from numerical analysis. 

The choice uf Uie size of ihe grid, which has to cover the region where plastic 
flow Is expected to occur, Is of utmost importance. A loo course grid will not detect 
changes In the values of plastic-strain for small loading increments. A too fine mesh 
size may result in distorted values of second-order derivatives of plastic strains, 
which appear in the function g(x,y). The loading Increment and Ihe grid size are re- 
lated to each other. A bud choice of either of them could result in the divergence of the 
iterative process. To allow tile maximum of grid points to be within the expected plus- 
tie zone, a variable grid spacing was chosen. The grid used for plane strain conditions 
wus finer, in general, than the one used for plane streSB cases. 

The interior region, where plastic flow is expected, was divided into r># rec- 
tangular cells. Due to symmetry alwut Ihe x-uxis, the number of unknown functions g, 
appearing in tlie bourn) y Eqs. (32) and stress Eqs, (36), Was reduced from r *• s to 
m - r > (s * l)/2, where now the coefficients of these functions represent the sum of 
the effect of left-hand and right-hand sides of the plastic Held. Because of computation 
tln.e limitations, (he grid was arranged in a 27 x 23 cell pattern, resulting in the num- 
ber of unknowns g to be equal to 324. By increasing the number of unknowns to 400, 


U’ 


the computation time for one iteration almost doubled. The smallest colls, located In 
the vicinity of the tip of the notch, have dimensions <5 /w * 0.004; 6 /w = 0.008 for 
plane strain cases, and $ x /w * 0.0C-1, 6 v /w = 0.016 for plane stress cases. 

The solution to the problem was obtained by tiro method of successive elastic solu- 
tions us discussed for the torsion problem and described in detail In Refs. 0 and (15) . 
Tlie computations wore performed on a digital computer using a FORTRAN IV program 
with single-precision arithmetic. The matrix system given by Eq. (35) was solved 
using the modified Gauss elimination method, which utilises pivoting and forward and 
backward substitutions. 


Results 


A number of beam problems were solved for both plane strain and plane stress 
cases. These Included notch depth to beam depth ratios of 0. 3 and 0. 5, notch angles of 
3° and 10 r , strain hardening parameter values of 0.05 and 0. 10. In addition, calcula- 
tions wore performed using the actual stress-strain curve of a 5083-0 aluminum alloy. 
For all cases Poisson's ratio was set at 0 

The loud increment size used was necessarily a compromise between the accuracy 
desired and computational time required for convergence. For strain hardening param- 
eter m = 0.05 the load increment size Aq defined as An max /o 0 was taken equal to 
0. 05; while for m = 0. 10, Aq = 0. 10. For the case of a 5083-0 aluminum alloy, where 
the actual stress-strain curve was used, the load was Incremented by Aq = 0. 025. 

For the beam with dimensionless notch depth a = 0. 5 the minimum load required 
to produce plastic flow at the most highly stressed grid points was found to be q = 0. 30, 
and for a * 0.3 the initial load was found (o be q » 0. So. The maximum load consid- 
ered was q = 0.7 for the a = 0.5 cases, and q = 0.9 for the a = 0.3 cases. In the 
process of solving the aforementioned problems, the case with strain hardening param- 
eter m = o. 0t> required approximately 50 iterations for each increment of load <1. o. , 

Aq = 0.05) lyr the relatively fine convergence parameter used. For cases where tne 
strain hardening parameter m = 0. 10 the average number of iterations needed for «ach 
increment of load (i. e. , Aq = 0. 10) was reduced to 40, while use of the actual stresa- 
strain curve resulted in convergence in approximately 10 iterations for the plane strain 
case and In 20 iterations for the plane stress case. 

Typical results of the compulations arc- presented In Figs, 11 to 18 and Tables V 
and VI. Complete detailed results are given in Ref. (11) . 

The growth of the plastic /one with load is shown in Figs. 11 to 14. It is seen that 
the shapes of the elasloplustie boundaries remain similar to each other us the load in- 
creases. As expected, plastic flow starts around the tip of the notch and ns the load 
increases appears also at the boundary opposite the notch. Comparison of Figs. 11 
and 12 with Figs. 13 and 14 shows that for the same loads the size of the plastic zones 
for plane strain arts considerably smaller than for plane stress. 

In the case of an clastoplastic problem the stress intensity factor Kj must be 
generalized t > the lorm 


KJttf 


max 



(37) 


where the exponent n is a function of the applied load, <r max . For linear clastic 


behavior Kf is identical with K { and n - 1/2. For the elasloplastlc case the varia- 
tion of n with load is shown in Tables V and VI. In the case of plane strain the stress 
singularity n decreases slowly as the load increases. For the plane stress case, 
there is a sudden drop in n from its elastic value as plastic flow appears. Subse- 
quently n slowly increases approaching a limit as tho load increases. 

Variation of the dimensionless generalized stress intensity factor with load is 
shown in Fig. 15 for the case of a specimen Hath notch depth of 5? = 0. G and a~ 10°, 
under plane strain condition and two values of strain hardening parameter m. The 
stress intensity factor shows no significant increase over the linear elastic value up to 
an applied load of q = 0.4l). Above this load KJ increases progressively for both 
m's, at tho faster rate for the lower strain hardening parameter, 

The products of y-di rocttonal stress and total strain were also calculated for var- 
ious eases. The order of singularity of that product was determined by platting 
Infer e ) against In r and by making a least squares fit of a straight line through the 
plotted points. It was found to be very close to unity for all cases considered. 

The y-dlrcctloual notch opening displacement was obtained for each case by numor- 
icrl integration of the relation = (l/2)(Uj ^ j) along straight line paths. For 
each ease a number of paths were chosen through the plastic region near the notch, and 
the resulting displacements were averaged. In general, tho notch opening displacement 
varies linearly with the load until the plastic zone is established at the boundary oppo- 
site the notch. Then it increases rapiaiy, reaching values several times that which 
would be calculated from the elastic solution. 

In order to verify in part the accuracy of the method used, a comparison of notch 
opening displacements was made with experimental results obtained by Bubsey and 
Jones (private communication from R. T. Bubsoy and M. H. Jones of NASA Lewis 
Research Center). The specimen used in this experiment, made of aluminum 5083-0 
with a length to width ratio of <1 and a crack length a = 0. 5, was subjected to three- 
point bending. The stress-strain curve for (ills specimen is shown in Fig. 10. The ex- 
perimental results as shown in Fig. 17 are in good agreement with numerical results 
obtained heroin. 

Finally, the J integral was evaluated for several cases. As in notch opening dis- 
placement calculations, straight lino paths wore chosen through the plastic zone near 
the tip of the notch. The integral was evaluated using values of stresses, strains, and 
displacements at cell centroids for a number of paths. The path independence of J was 
not conclusive, since tho results varied up to 15 percent from the averaged value. It is 
possible that the results obtained herein do not indicate that the path independent prop- 
erty is lost but rather that the field values of tile displacements are not calculated with 
sufficient accuracy. 

The average values of the dimensionless .T integral as a function of load arc 
plotted in Fig. 18 for a ease of a specimen with u 10° edge notch, a = 0.5, m — 0.05, 
and plane strain condition. A! the start of plastic flow' J increases rapidly with load. 
This is followed by almost linear nation with additional load. 

From the above results it , * oars that the BIE method applied to the plane problem 
and formulated in terms of the \l.j stress function is capable of giving detailed results 
such as stress and strain distributions around the tip of die notch and, related to them, 
the shapes of plastic zones. This was accomplished using a relatively small number of 
unknowns . 
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The presence of a singularity at the Up of the notch makes accurate answers very 
difficult to obtain. Nevertheless good agreement was obtained between the calculated 
results and ejqjerimcn tally measured notch opening displacement os shown in Fig. 17. 
Some improvement in the solution techniques and further invesUgaUon of die Influence 
of the boundary nodal spacing and interior grid slxe on the resulting stress and strain 
fields, and therefore, on the notch opening displacements and J integrals, is still 
desirable. 

THE DISPLACEMENT FORMULATION 


Although the biharmonic formulation previously described appears as a viable 
approach in solving the plane elastoplastlc problem, some difficulties are encountered 
in calculating displacements since numerical differentiations are required in the proo- 
cess which can lead to appreciable error* and inconsistency of results, A more direct 
formulation of the problem is given in terms of the Wavier equilibrium equations for the 
displacements. The general equations are given in Ref. (4) and several problems using 
these relations are reported in Ref. (16). 

As shown in Ref. 0 the Navier equations with plastic flow can be converted to the 

B1E 

AU i( P) (UjjPj - T IjUj )dq + S JW (rj> k ♦ *f k )dA (38) 

where Uj and Pj are the boundary displacements and boundary loads, respectively, 
and the usual tensor notation is used. The tensors U^, T^, and are given by 



U ij “ C l (6 ij C 2 ln r " r * i r » j> 

(39) 

T - 2 * 

T ij~7 

(^i) C 4 + 2r , i r , )) + C 4 (r, jn t - r,^) 

(40) 


_q J 





V - ± 

w Jhi~ r 

^4 <6 ij r 'k + 6 ki r * j ~ 6 Jk r » + 2r ‘ i r * j r *k] 

(41) 


with 


C l = 


1 _ 

8jtG( 1 - p)' 


-, C« = 3 - 4p 


: 3 * — . C 4 = 1 - 2p 


(42) 


and r is the distance from the fixed point P to the variable point of integration, q. 
The above equations are for the case of plane strain. For plane stress one replaces 
by p/(l + p) . The coefficient \ is equal to 1, if P is an interior point and is equal 
1/2, if P » p is a boundary point. 

The solution is now obtained by replacing the integrals by sumB as before, result^ 
ing in 2n simultaneous equations for the case of n boundary segments. These can be 
written as the matrix equation 


S' ■*= 


A+-M B u y C AC 

2 

, = + + (43) 

a* B * f - 1 lyj L a J L d J Ia d . 

whore A, B, A', and B' arc n «• n matrices of known coefficients. u and v are the 
x and y displacement vectors for the boundary nodal points. > and 6 are vectors 
given by 

y = oP x + iP y j 

? (44) 

6 ~ a'P x * 0*!» 

J 

where &, &, a', and d' arc known n ■ n matrices and P x , P y are the x and y 
boundary force vectors at the nodal points. The vector# C and D are functions of 
the accumulated plastic flow and are known at the beginning of each load increment. 

The vectors AC und AD dopond on the plastic flaw increments during the current 
load increment and are obtained by iteration. All the terms appearing in Eqs. (43) 
and (44) are listed in Appendix C. 

For the first boundary value prohlem where the loads are specified over the com- 
plete boundary, Eq. (13) is solved directly for the unknown displacements at the bound- 
ary nodes. For a mixed boundary value problem, where the displacements at some of 
the nodes arc specified and the loads are unknown at these nodes an obvious interchange 
of the appropriate columns of the coefficient matrices must be made. 

Once Eq. (43) is solved for the unknown displacements (and loads if any), the total 
strains are computed at any interior point by differentiating Eq. (38) with a = 1. At an 
interior point Pjj s P(Xj,y^) we can write 


t x (P ij > ~ ^ (P xm E ijm + P ym 1- ijm “ u m G ijm “ v m I, ijm > + *ij + AI ij 
m=l 

n 

f y (P ij ) = X) (P xm E ljm 1 P ym F Ijm “ u m G ijm " v m H ijm* + *lj + 


y< p ij> - S (P xm E ij * p ym F ijm “ u m G ijm ' v m»ljm> + <ij + *\\ 


All the coefficients arc given in Appendix C. The plastic strain increments are then 
computed from the plastic strain- total strain relations and the stress-strain curve ns 
given in Refs. (3,j5), the method of successive clastic solutions being used to obtain 
these strain increments iteratively. 

It is shown in Ref. <1(>), that greater accuracy can be obtained for the same number 
of nodal points by assuming linear variations of the unknowns on the boundary Intervals. 
This of course complicates to some extent the calculation of (he coefficient matrices. 
Furthermore the nodal points can no longer bo taken at the midpoints of the intervals, 
but must be taken at the end points. This introduces some further complication when a 
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nodal point occur* at a comer, since firstly, the boundary load may be discontinuous at 
a comer and secondly, the jump in the boundary integral at a corner is no longer 1/2, 
but depends on the comer angle. 

Although these additional complications can be taken care of as was done in 
Ref. (Jj6), a somewhat different but approximately equivalent approach was attempted 
herein. The nodal points were kept at the centers of the Intervals. Each Interval was 
divided into a number of subintervais and the integral for each subinterval was weighted 
linearly according to its distance from the nodal point. This approximates a linear var- 
iation of the unknowns over the intervals. 

This approach did indeed give improved results. For example, forthe beam of 
Fig. 7, with a/w = 0.5, L/w = 1.2, the plane strain elastic maximum crack opening is 
Ev/6M = 5.67. The value obtained using 60 nodal points was 4.92 assuming constant 
values over each interval, and was 5. 62 using the linear weighting technique. 

Results 

Calculations were made by the above technique for the same problem as was solved 
by the biharmontc formulation previously described. Some preliminary results are 
shown In Figs. 17 and 19. Fig. 17 compares both the biharmonic fomulation and the 
displacement formulation with experimental results for the maximum notch opening. 

The beharmonic and displacement formulations give results which are in very good 
agreement. The same is seen in Fig. 19 where the stress at a value of x very 
close to the notch tip is plotted against the distance from the notch centerline, v. The 
agreement between the two formulations is again very good. The calculations for larger 
load Increments using the displacement formulation have not yet been carried out. 

These are presently under way. 

A preliminary comparison of the convergence rate of the two formulations indicates 
(hat the plastic flow computations converge more rapidly using the displacement formu- 
lation. A comparison of the computer times, however, could not be made, since the 
two types of calculation were carried out on different computers. 

CONCLUDING REMARKS 

This preliminary Burvey of the use of BiE methods for elastoplastlc problems indi- 
cates that they form a viable und worthwhile approach for solving such problems. The 
torsion problem in particular can easily be solved for almost any geometry cross 
section. 

The plane elastoplastic problem can be solved by using either a blharmonic formu- 
lation or a displacement formulation. Both appear to give good results with relatively 
small sets of equations, even for problems with singularities, such as beams with 
notches. 

Althou^i no comparison was made herein with the finite element method, such com- 
parisons were made in Ref. (16). It is indicated in Ref. (16) that the computer times 
for the finite element method and BIE method (using the displacement formulation) are 
comparable. Finer resolution can however be obtained by the BIE method. 

The application of the BiE method to elastoplastic problems is still in its early sta- 
ges. Much work remains to be done in refining the techniques for optimum application. 
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APPENDIX A 


BOUNDARY INTEGRAL COEFFICIENTS FOR TORSION PROBLEM 

The division of the boundary into intervals with their corresponding nodal points is 
shown in Fig. 30. The x and y coordinates of a boundary nodal point pj are desig- 
nated as (Xjjj, y^) . Hie coordinates at the beginning and end of an interval (say inter- 
val j) are designated by (i j, tij) at the beginning of the interval and by ((jj +1 , tj« + 1 ) at 
the end of the interval . The interval lengths lij need not be equal. The coordinates of 
the centroid of on interior cell where plastic flow occurs are designated by (x^, y ( ). 

The coofflcionts in Eq. (13) are then given by 


W/2) 


F(l J-(l/2) 


— (to r a )*l 

% W 


<*bj " x bl> £ j * <>bi - >'bi» r, 'j 


“u " v ” A a ik 


The last relation follows from the Gaussian condition, that is. 


f ~ In 

l/C <i 


r pq dqa ” r 


To evaluate the b^ coefficients Simpson's rule is used for the case i 4 j and 
closed form integration is used for tile case i = j since the integrand is singular for 
1 *> j. The result is 

b l|*^' n [ r i.J-<l/2) V'" r UHl/2)J 


•*11 • h i (’" 7 - ') 


R ‘~^ 


kf ln r ikf ^kf 


where ^ is the sum for all the plastic cells in the region and AA^ is the area of 

M 

die coll with coordinates (x^, y f ). 

The coefficients Ay^, C^, and Djj kf are given as follows using Simpson's 


rule: 





APPENDIX D 


COEFFICIENTS OF THE 8TREB8 EQUATIONS (36) 

The coefficients appearing in stress Eqs. (36) are given by the following relations: 



(continued) 
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APPENDIX C 


COEFFICIENTS OF DISPLACEMENT FORMULATION EQUATIONS 
The coefficient* appearing in Eqs. (43) to (45) are given aa follow* 


Let 


whore 




m+1 

m 


E g = sin 0 co* o 


m+1 


E 3 * ain‘0 


E 5 - In r| 


e 6 - «ns*e] 


m 

m+1 

m 

m+1 

Im 

m+1 

m 


E, - *in*0 


m+1 

m 


E - tan 0 In r| 

o 


m+1 

m 


E g “ tan 0 


m+1 

m 


D = r cos U 
m m 


0_ = sin' 
m 


■*(I«rsLi2 

V r im 


m *i m x m ” x 
” m m— ' 
r im 


‘) 


(54) 


Then 


Ajj a C 3 [(C 4 + l)Ej + (l 2 “ - SlmEijj 

B tJ - C s [2lmE 2 + (I 2 - m 2 )E 3 + c 4 eJ 

A|j - CgfaimE, + <t 2 - m 2 )^ - C 4 E 8 ] 

% * c s[(C 4 + 1)^! - («' - m^Eg +SlmE^ 

0 ij " Dc l[ C 2 " < <2 ‘ - <C 2 + m^Eg + 2lmE 5 + CgE g J 

hi " °Ij " DC 1 [- 24 «8 1 + lm * 9 - <** - <n*>E 5 ] 

= CjD[(C 2 + i 2 - m^Ej - (C 2 + (^Eg - 21m Eg + CjE^ 

E ljm “ ' c l[ (C 2 ” (£2 " m2 »* E i + ( f3 - 31m 2 >E 2 

+ (m 3 - 3ml 2 )E 3 . m(Cg - Sl^Eg] 

F i]m = " Cj^^mEj - (m 3 - 3t 2 m)E 2 + (l 3 - Slm^Eg - 1(1 2 - m^EgJ 
E ijm = - C l[“ 

F ljm = - C l[< 

- (I 3 - 3tm 2 )E 3 + 1(C 2 - 2m 2 ) E 5 j 

E i]m = “ C l[ <C 4 ‘ ({2 - m2 » mE ! - ( m ' J - 31 2 m)E 2 

+ (l 3 - 31m 2 ) Eg + 1(C 4 ♦ 2m 2 )Eg| 

F ijm “ ’ c l[ lC 4 " f2 ~ m2 > CE l " < f3 * : >«m 2 )E 2 

+ (m 3 - 31m 2 )E 3 - m(C 4 + 2t 2 )E 5 J 

£ 

G ijm = — ~ [<• " 6 *“ " c 4 > f % + (C 4 “ 2 + 10{2)rnE 3 

+ 21(41 2 - $)E 4 + 2m (4m 2 - 3 )Et| 


21 m 2 Ej - (l 3 - 31m 2 ) E 2 - (m 3 - 3l 2 m)Eg - m(l 2 - m 2 )E 8 ] 


(C 2 + l 2 - m 2 )mEj + (m 3 - 31 2 m)E 2 


(5S) 


> (06) 


(continued) 


»** 


> (M) 


0 

H ljm 6<8 “ C 4> mE 2 + <" C 4 • 3 * 1 °w»Vb 3 

- 2m(4m 2 - 3)B 4 + 2£(4£ 2 - 3)E,J 

GJj m = J(2 - 6m 2 - C ( )£E ? - (-C, - 2 * 10( 2 )inK 3 

- 2£(4£ 2 - 3)E 4 +■ 2m<4m 2 - 3)E ? j 

Hljm “ -^[< 4 ’ flm2 * C 4> mE 2 - (C l - 2 4 l° m2 > £ E 3 

+ 2m(4m 2 - 3)£ 4 - 2l(4t 2 - 3)E^J 

Glj m a jm(6m 2 - 5)E 3 + I(l0m 2 - 3)E 3 - 2m(4m 2 - 3)E 4 + 2£(4£ 2 - 3)E-j 
|||j m «-^jf ( 6£ 2 - 5)E 2 + n>(3 - 10( 2 )E 3 - 2£(4£ 2 - 3)E 4 - 2m(4m 2 - 3)E 7 J 


and 

‘ij 31 [ 8 lUl (P iJ ,Q k£ )t x <Q k£ ) + 8 ll22 (p lj- Q k£> t £ (C W 


+ 2 s U12 (P tJ* Q k£ )t xy (Q k£ ) ] AA W I 


l ij 3 £ kt (^l^x f S 2222'y 4 2 '^l^xyj^kt 
l lj “ S k( ( S 121l c x 4 S 1222 e y 4 2 S 12l2*xy) AA k£ 


(57) 


und 



(53) 


where k£ represents an interior cell. 

If the nodal point i and the Integration Interval j are on the same straight seg- 
ment of the boundary, then D - 0 and </ - iir/2, where the plus sign Is used when j Is 
ahead of i, and the minus sljm Is used when j Is behind 1, This leads lo 0** In the 
calculation of dy and dfy. A simple limiting process shows that for this case, with 
i ^ j. 
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TABLE J. - COMPARISON OF VALUES Or DIMENSIONLESS WARPING FUNCTION ON BOUNDARY 


OF ELASTIC SQUARE PLATE WITH EXACT ELASTIC SOLUTION 


Boundary 

value, 

y 

Exact 

warping 

function, 

w 

Value of warping function by 
boundary integral method 

Boundary 

value, 

y 

Exact 

warping 

function, 

w 

Value of warping function by 
boundary intogral method 

Interval*, n 


Intervale, n 


m 

m 

12 

16 

4 

8 

12 

16 

0.03125 

0. 01095 




0.01095 


0. 1424 




0. 1424 

. 04167 

.01459 



0.01459 


.5417 

. 1433 



0. 1433 


.0525 

.02185 


0.02184 



.5625 

.1448 


0. 1446 



.00375 

. 03264 




.03264 

.5938 

.1461 




. 1460 

. 1250 

. 04328 

0. 04311 


. 04328 


.6250 

.1461 

0. 1446 


. 1461 


. 1563 

.053T4 




.05372 

.6563 

.1449 




. 1449 

.1875 

.06380 


.06386 



.6875 

.1422 


.1419 



.2083 

.07060 



.07051 


.7083 

. 1396 



.1395 


.2188 

. 07380 




.07377 

.7188 

.1380 





2813 

.09235 




.09234 

.7813 

.1244 




iH 7TV 

2017 

.09527 



. 09525 


.7917 

.1214 



.1212 


.3125 

.1009 


.1009 



.8125 

.1148 


. 1140 



. 3838 

.1090 




.1090 

.8438 

.1029 




.1027 

.3750 


.1159 


.1164 


.8750 

. 06864 

.08826 


.08811 


. 4063 

.1233 




.1232 

.9063 

.07166 




.07129 

.4375 

B 


. 1282 



.9375 

.06160 


.06142 



.4583 

1 



.1329 


.9583 

.03644 



.03621 


.4688 

. 1346 





. 1346 

.9688 

. u.1808 




.02706 
































































































TABLE V. - ORDER OF STRESS SINGULARITY n AT the TIP OF T1IE NOTCH 
FOR A SPECIMEN WITH A SINGLE EDGE NOTCH SUBJECTED TO 
PURE BENDING • PLANE STRAIN 


[Poisson's ratio u 0.33] 


Dimension- 
less notch 
depth, 
a 

Notch 

angle. 

*r. 

top 

Strum Imrcl* 

1‘lUlHi |>4* 

m 

Elastic 


Dim 

onsionlcss load. q 


0 4 

0.5 

0.0 

■ 

OH 

0.9 

0.3 

3 

■K 

0.40MB 

Mj 

0.41(0 

0.490 

0,487 

0 473 

0.475 

.3 

10 





490 

497 

492 

.480 

487 

6 

10 




0 490 

. 400 

480 

472 


BBJH 

5 

10 

MUM 




406 

400 

400 

.478 

— 

BBS 


TABLE VI. - ORDER OF STRESS SINGULARITY n 
AT THE TIP OF THE NOTCH FOR A SPECIMEN 
WITH A SINGLE EDGE NOTCH SUBJECTED 
TO PURE BENDING - PLANE STRESS 


[Dimensionless notch depth a - 0.3; notch angle 
i» lo' J , strain hardening iiaramctor 111 -0.10; 
Poisson’s ratio tj 0. 33. J 


Ll.iHUr 

Dimensionless load, q 

0 5 

0 6 

0.7 

0. 0 

0. 0 

0 4999 

0.419 

0 434 

0. 448 

0.451 

0. 450 
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Figure?. -Region R, boundary curve c, and 
geometric quantities entering Into boundary 
integrals. 




Figure 4. - Variation of dimensionless moment with dimensionless 
angle of twist per unit length for several values of strain- hardening 
parameter for square cross section. 
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Rgure 5. - Variation of dimens ionless maximum shear stress with 
dimensionless angle of twist per unit length for several values of 
strain- hardening parameter for sguare cross section. 
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Figure 11. - Growth of plastic zone size with load for specimen with single edge notch subjected to pure 
bending, Plane strain; dimensionless notch depth a * 0. 3; notch angle c ■ 10°; strain hardening 
parameter m -0.10; Poisson's ratio p-0.33. 
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Figure 1?. - Growth of plastic zone size with load in vicinity of notch for specimen 
with single edge notch subjected to pure bending. Plane strain; dimensionless 
notch depth a *0.3; notch angle a* 10°; strain hardening parameter 
m • 0. 10; Poisson's ratio p • 0. 33. 





Figure 13. - Growth of plastic zone size with load for specimen with single edge notch subjected to pure 
bending. Plane stress; dimensionless notch depth I ■ 0. 3; notch angle a ■ 10° ; strain hardening 
parameter m ■ 0. 10; Poisson's ratio \i • 0. 33. 
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Figure 14. - Growth of plastic zone size with load in vicinity of notch for a specimen 
with single edge notch subjected to pure bending. Plane stress; dimensionless 
notch depth a * 0. 3; notch angle a • IQP; strain hardening parameter m • 0. 10; 
Poisson's ratio p « 0.33. 
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Figure 1$. -Variation of dimensionless generalized stress intensity factor with load for 
specimen with single edge notch subjected to pure bending. Plane strain; dimensionless 
notch depth t * 0. 5; notch angle a • 10°j Poisson's ratio |i * 0. 33. 
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Figure 16. - Stress-strain curve for aluminum 5083-0 used in test (private communication from 
R. T. Bubsey and M. H. Jones of NASA Lewis Research Center). Young's modulus of elas- 
ticity E ■ 7. 79x10 s newtons per square centimeter (11. 3x10$ lb/in. 2 ); Poisson's ratio 
U-0.33. 
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Figure 18. - Dimensionless plane strain J Integral for specimen with single edge notch 
^ subjected to pure bending. Dimensionless notch depth a • 0. 5j notch angle o • 10°; 

strain hardening parameter m * 0. 05; Poisson's ratio y ■ 0. 33. 
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Figure 19. - Dimensionless y- 
directlonal stress at a function of 
distance from notch centerline, 
x • D 0074, a • Ol 5, q * 0. 4. 
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Figure 20. - Boundary and interior notation for computing «■ 
efficients given in Appendix A. 
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